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NOMENCLATURE 


wall  friction  coefficient  «  - y 

1  *  * 


wall  pressure  coefficient 


1  *  * 

JP 


spline  derivative  approximation  of  ^ 


spline  derivative  approximation  of  ■*— 


spline  derivative  approximation  of  — y 

9n 


spline  derivative  approximation  of  — j 

3n 


reference  length  (dimensional) 


1 
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transformed  normal  coordinate 

step  size  in  n-direction 

wall  static  pressure  (dimensional) 


Reynolds  number  ■  — j— 

V 

transformed  axial  coordinate 
step  size  in  s-directioo 


free-stream  speed  (dimensional) 


-6- 


25  January  1982 
GHHrmmj 


axial  coordinate 

axial  coordinate  of  initial  line  in  computational  domain 
axial  coordinate  of  end  of  region  I  in  computational  domain 
axial  coordinate  of  end  of  region  II  in  computational  domain 
axial  coordinate  of  end  of  computational  domain 
normal  coordinate 

normal  coordinate  of  outer  edge  of  computational  domain 
normal  coordinate  of  wall 

first  derivative  of  wall  coordinate  with  respect  to  x 
second  derivative  of  wall  coordinate  with  respect  to  x 
normal  coordinate  of  front  plate 


normal  coordinate  of  rear  plate 
velocity  component  in  x-direction 
velocity  component  in  y-direction 


boundary-layer  thickness  at  x  =  x 


0 


stream  function 


vorticity  magnitude 

step  size  ratio  in  n  direction — see  Eq.  (42) 
Blasius  similarity  variable — see  Eq.  (100) 
wall  shear  stress  (dimensional) 


kinematic  viscosity  (dimensional) 


fluid  density  (dimensional) 


All  other  quantities  are  defined  in  the  text. 
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All  quantities  in  the  text  are  made  dimensionless  as  follows: 

* 

distance  by  L 

* 

velocities  by  U_ 

OO 

it  it 

stream  function  by  L  U 

it  it 

vorticity  by  U^/L 
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I .  INTRODUCTION 


The  motivation  for  the  present  study  is  the  development  of  an  efficient 
numerical  calculation  procedure  for  flows  at  moderate  to  high  Reynolds  numbers 
where  a  strong  interaction  exists  between  a  boundary  layer  and  a  surrounding 
vorticity-free  flow.  One  possibility  for  achieving  efficiency  might  be  to 
solve  a  set  of  equations  which  is  simpler  than  the  Navier-Stokes  equations. 

In  this  regard,  a  number  of  papers  have  been  published  in  recent  years  [5-8]* 
which  show  that  strong  interaction  flow  fields  may  be  calculated  with 
considerable  accuracy  by  using  an  approximate  set  of  equations  referred  to  as 
the  parabolized  or  thin-layer  Navier-Stokes  equations. 

Four  parabolic  models  of  the  Navier-Stokes  equations  are  possible, 
according  to  Davis  and  Rubin  [1].  At  the  top  of  this  hierarchy,  and  the 
one  addressed  in  this  report,  is  the  parabolized  vorticity  equation  approxi¬ 
mation.  In  this  approximation,  the  Navier-Stokes  equations  are  written, 
preferably,  in  an  optimal  or  nearly  optimal  coordinate  system.  Then  all  terms 
that  make  the  vorticity  equation  elliptic  are  neglected,  while  all  terms  in  the 
elliptic  stream  function  equation  are  retained  in  order  to  allow  for  a  proper 
interaction  between  the  viscid  and  inviscid  layers. 

As  discussed  by  Davis  [2],  the  accuracy  of  a  parabolic  model  depends  on 
the  choice  of  coordinates.  Thus,  optimal  coordinates  in  the  sense  of  Kaplun 
[3],  provides  the  best  accuracy  for  the  parabolic  model,  as  demonstrated  in 
the  calculations  of  Davis  and  Rout  [4].  The  determination  of  optimal 
coordinates  is  itself  a  formidable  task  because  a  complicated  set  of  coupled 
flow  and  coordinate  equations  must  be  solved.  Studies  of  the  parabolized 
vorticity  model  using  conformal  coordinates  (which  are  nearly  optimal)  have 
been  made  by  Davis  [5],  Ghia  and  Davis  [6]  and  Napolitano,  Werle  and  Davis 
[7].  These  results,  even  when  separated  zones  are  present,  are  nearly 
identical  to  full  Navier-Stokes  solutions  at  intermediate  and  high  Reynolds 
numbers.  The  work  of  Steger  [8]  for  high  Reynolds  number  compressible 
viscous/ turbulent  flow  about  an  airfoil  using  the  parabolized  vorticity 
approximation  indicates  that  accurate  results  can  even  be  obtained  for 
numerically  generated  grids  that  are  slightly  non-orthogonal  and  not  chosen 
from  optimality  considerations. 

The  present  work  provides  a  further  test  of  the  accuracy  of  the  parabolized 
vorticity  approximations  in  a  non-optimal  coordinate  system  as  applied  to  the 
calculation  of  strongly  interacting  flows  at  intermediate  to  high  Reynolds 
numbers.  Toward  this  end,  the  following  problem  is  solved:  A  steady, 
incompressible,  two-dimensional,  laminar  flow  exists  over  two  flat  plates 
aligned  with  the  oncoming  stream,  but  separated  vertically  and  horizontally. 

The  plates  are  joined  by  a  smooth  center  section  described  by  a  cubic  so  that 
the  wall  has  continuous  slopes  at  the  junctures.  The  pressure  distribution 
in  the  boundary  layer  thus  may  be  varied  by  changing  the  vertical  and  horizontal 
spacing  of  the  fore  and  aft  plates. 


Numbers  in  brackets  []  indicate  References,  see  Page  41. 
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The  computational  domain  for  this  problem  is  shown  schematically  in 
Fig.  1.  On  the  initial  surface  at  x  =  x  ,  measured  from  the  leading  edge 
of  the  first  plate,  a  Blasius  boundary  layer  exists  between  y  =  0  and  y  =  6 
(zone  1).  Between  y  =  <SQ  and  y  =  y  ,  the  outer  edge  of  the  computational  ° 
domain,  a  uniform  flow  is  present  (zone  2).  To  ensure  that  a  strong  inter¬ 
action  takes  place  between  the  boundary  layer  and  the  outer  inviscid  flow, 
zone  1  is  taken  to  be  a  significant  fraction  of  zone  2  and  the  line  y  =  y 
is  taken  to  be  a  streamline.  Actually,  the  problem  posed  is  the  lower  e 
half  of  a  diffuser  flow  that  is  not  fully  developed. 


The  method  used  here  to  generate  a  non-optimal  wall-fitted  coordinate 
system  is  to  apply  a  shearing  transformation  to  the  computational  domain, 
shown  in  Fig.  1,  which  maps  this  domain  into  a  rectangle.  The  vertical 


distance  between  the  two  plates,  y, 


W1 


-  y, 


W2 


is  restricted  to  be  much  less 


than  the  horizontal  distance,  X2  -  x^,  to  keep  separated  flow  zones  thin 
so  as  not  to  violate  the  parabolic  approximation.  Consequently,  the 
coordinate  system  is  only  mildly  non-orthogonal .  A  typical  sheared  grid 
is  shown  in  Fig.  2. 


The  numerical  solution  strategy  adopted  to  solve  the  above  problem 
is  along  the  lines  of  Murphy  [9]  in  his  development  of  an  efficient 
Navier-Stokes  solver.  The  main  features  are: 


1.  The  steady  parabolized  vorticity  model  is  solved  using  line 
overrelaxation  in  conjunction  with  the  Newton-Raphson  technique. 

2.  Fourth-order  spline  discretization  is  used  in  the  direction  normal 
to  the  main  flow  so  that  the  region  of  high  flow  gradients  near 
the  wall  can  be  resolved  with  relatively  few  grid  points.  Standard 
second-order  finite  difference  discretization  is  then  used  in  the 
main  flow  direction  where  flow  gradients  are  expected  to  be  mild. 

The  classical  Peaceman-Rachford  ADI  procedure  [10],  used  in  References 
4-7,  requires  an  orthogonal  coordinate  system  for  stability.  For  non- 
orthogonal  coordinates,  a  more  general  splitting  technique,  such  as  that 
of  Beam  and  Warming  [11],  must  be  used.  On  the  other  hand,  SLOR,  the 
procedure  used  here,  is  a  straightforward  and  well-proven  alternative  for 
obtaining  steady  solutions  which  is  not  affected  by  mild  non-orthogonality 
of  the  coordinate  system. 

Spline  collocation  techniques  as  a  means  of  discretizing  viscous  flow 
problems  have  been  extensively  studied  by  Rubin  and  Khosla  [12,  13].  They 
point  out  the  following  advantages  of  the  spline  formulation  over  three  or 
five  point  finite  difference  discretization: 
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1.  The  discretized  equations  formed  with  splines  remain  block 
tridiagonal . 

2.  The  treatment  of  derivative  boundary  conditions  is  more 
straightforward  because  the  derivatives  are  treated  as 
unknowns . 

3.  A  highly  nonuniform  grid  produces  less  deterioration  in 
accuracy  than  with  finite  differences. 

4.  For  equal  accuracy  the  fourth-order  spline  requires  one- 
fourth  the  number  of  grid  points  in  a  given  direction 
compared  to  second-order  finite  differences. 

The  disadvantages  of  splines  are  that  the  discretized  equations 
are  more  complicated  than  their  finite  difference  counterparts  and 
the  storage  requirements  are  larger  for  a  given  number  of  mesh  points. 
Since  fewer  points  are  required  for  comparable  accuracy,  the  spline 
approach  is  still  the  more  efficient. 

The  present  discretization  is  similar  to  that  of  Murphy  [9],  who 
used  a  generalized  Galerkin  technique  with  splined  Taylor  series 
expansions  to  achieve  fourth-order  accuracy  in  the  direction  normal 
to  the  main  flow.  The  use  of  fourth-order  splines  achieves  essentially 
the  same  result  but  with  less  algebra  to  arrive  at  the  numerical 
algorithm. 

This  report  covers  numerical  solutions  using  the  parabolized 
vorticity  approximation  as  well  as  boundary-layer  solutions  to  show 
the  magnitude  of  the  strong  interaction  effect.  A  later  report  will 
present  comparisions  of  full  Navier-Stokes  solutions  with  the  parabolized 
vorticity  model. 

II.  ANALYSIS 

Governing  Differential  Equations 

The  starting  point  of  the  analysis  is  the  steady,  incompressible 
Navier-Stokes  equations  written  in  Cartesian  coordinates.  In  dimension¬ 
less  stream  function-vorticity  form  these  equations  are: 


Stream  Function 


Ip  +  if)  =  -  c, 
xx  yy 


Vorticity 


*y  ?x  "  ^y’i  Uxx  +  V 


(1) 


(2) 
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where  subscripts  denote  partial  differentiation  with  respect  to  the  subscripted 
variable.  For  the  truncated  computational  domain,  shown  in  Fig.  1,  the  boundary 
conditions  are  similar  to  those  of  Murphy  [9].  At  x  =  x^,  the  flow  is  pre¬ 
scribed  : 


ip  =  <Kx0,y)  ,  (3a) 

C  =  ;(xQ,y)  .  (3b) 

At  y  =  yw(x)  the  no-slip  and  impervious  wall  conditions  prevail: 

u  =  v  =  0  at  y  =  yw(x) 
which  are  equivalent  to 

<Kx,y  )  =  0  ,  (4a) 

w 

ty(x,yw)  =  0  .  (4b) 

At  y  =  ye,  the  outer  edge  of  the  domain,  the  stream  function,  and  vorticity 
are  prescribed: 

iMx,ye)  =  ^e(x)  *  (5a> 

C(x,ye)  =  Cg(x)  ,  (5b) 

where  usually  £e  =  0.  At  the  downstream  boundary,  x  =  x^,  the  elliptic  region 
is  terminated.  Thus 

*xx(*f.y)  =  0  >  (6a) 

5xx(xf,y)  “  0  *  (6b) 


These  last  two  conditions  are  equivalent  to  requiring  the  flow  to  satisfy 
boundary-layer  equations  at  the  outflow  boundary. 

To  produce  a  wall-fitted  coordinate  system,  the  following  shearing 
transformation  is  used: 

s  =>  |  [1  +  y^2]  1/>2  dx  ,  (7) 


n  *  $(x)  •  y  -  A(x) 


(8) 
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where : 


4>(x)  = 


ye  -  yw(x) 


yw(x) 

A(x)  =  - 7— r-  =  <P  y  .  i 

ye  "  yw(x)  w 

and  the  prime  in  Eq.  (7)  denotes  differentiation  with  respect  to  x.  Thus, 
s  corresponds  to  arc  length  along  the  wall  and  n,  the  normal  variable, 
ranges  between  0  and  1  as  y  ranges  between  yw  and  yg.  The  computational 
domain  is,  therefore,  mapped  into  the  rectangle  s^  <  s  _<  sj,  0  <  n£l. 

The  transformed  Navier-Stokes  equations  in  (s,n)  variables  are, 
noting  that  s  =  s  =0: 

y  yy 

Stream  Function 

2  2  2 

s  +2sniJ>  +(n  +n)^  +s  .  U  .  .  _  n 

x  ss  x  x  sn  x  y  nn  xx  s  +  ii  +  C  =  0 

J  xxTn 


Vorticit} 


s  n  (iji  C  -  C  )  *  V-  [s  +  2s  n  5  +  (n  2  +  n  2)  X, 

x  y  yn  s  s  n'  Re  x  ss  x  x  sn  x  y  nn 


+  s  Q  +  n  C  ]  •  (12) 

xx  s  xx  n 

The  metric  coefficients  sx,  nx,  etc.,  obtained  from  Eqs.  (7)  and  (8),  are 
derived  in  the  Appendix. 

The  parabolized  form  of  the  vorticity  equation  is  obtained  from  Eq.  (12) 
by  dropping  £gs  and  Sgn.  To  be  consistent,  the  term  sxx  £s  is  also  dropped. 
The  result  is: 

Parabolized  Vorticity 


8  n  (f  ^  -  ty  Z  )  -  ir~  [(n  2  +  n  ^)c  +  n  C  ] 

x  yvr  Ss  Ys*n  Re  x  T  y  snn  xx^nJ 


(13) 
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The  flow  in  the  computational  domain,  s„  s  s,,  0  <_  n  <  1,  is  governed 

by  Eqs.  (11)  and  (13). 

The  boundary  conditions  in  (s,n)  coordinates  are  as  follows: 


4/  =  iKs0,n)  ,  (14a) 

C  =  S(s0,n)  .  (14b) 

on  the  wall,  n  =  0: 

♦(s.O)  -  0  ,  (15a) 

Vs ’0)  =  0  ‘  (15b) 

On  the  outer  edge,  n  =  1: 

♦  (s ,  1)  *  <J>e(s)  ,  (16a) 

C(s,l)  =  5e(s)  (16b) 


On  the  downstream  boundary,  s  =  s^,  the  vorticity  equation  requires  no 
boundary  condition  because  it  is  parabolic.  The  only  boundary  condition 
required  is 


4»ss(Sf,  n)  -  0  , 

which  applies  to  the  stream  function  equation. 


(17) 


Spline-Finite  Difference  Equations  in  the  Field 


To  use  spline  approximations  in  the  n-dlrection,  the  following  spline 
derivatives  must  be  defined: 


£ 

L 


♦ 

♦ 


£ 


5 


L 


C 


nn 


nn 


(18) 

(19) 

(20) 
(21) 


The  stream  function  and  parabolised  vorticity  can  then  be  viewed  as  equations 
from  which  to  determine  1$  and  L^.  From  Eq.  (11)  one  obtains 
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and  from  Eq.  (13), 


where : 


L5=  _d  +  f (iL^4  -  JlS  ) 
s  s 


(23) 


c  = 


d  = 


2  s 

n 

X 

X 

2  J 

2 

> 

n  + 

n 

X 

y 

2 

s 

X 

2  , 

2 

> 

n  -f 

n 

X 

y 

s 

XX 

2 

2 

i 

+ 

e 

n 

X 

y 

a 

XX 

2  ^ 

2 

> 

n  + 

n 

X 

y 

1 

2  ^ 

2 

> 

n  + 

n 

X 

y 

Re  s 

n 

> 

:  y 

2  , 

2 

n  + 

n 

X 

y 

(24) 


(25) 


(26) 


(27) 


(28) 


(29) 


Note  that  the  cross  derivative  >|;sn  in  Eq.  (11)  has  been  treated  as  £ 


Equations  (22)  and  (23)  are  now  discretized  in  the  s-direction  as 
follows: 

1.  A  constant  step  size  As  is  used  for  simplicity. 

2.  In  the  stream  function  equation  (L^)  ,  second-order  accurate  central 
differences  are  used. 


3.  In  the  parabolized  vorticity  equation  (I/*),  three-point  backward 
differences  are  used  to  maintain  stability  (9], 
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At  line  i,  the  appropriate  finite  difference  expressions  for  s-derivatives 
to  be  used  in  Eq.  (22)  are 


<Vi 


i+1 


-  F 


2  As 


—  +  0(As2) 


(30) 


ss  i 

and  in  Eq-  (23), 


F.  ,  -  2F.  +  F  ,  „ 

(F__)  j  =  — - - - —  +  0(As2) 


As 


(31) 


F.  -  4F.  .  +  3F.  0 

(Vi -  - - +  0(As > 


(32) 


Considering  quantities  at  node  point  (i,j),  corresponding  to  (s.,n.) 
as  unknowns,  Eqs.  (22)  and  (23),  with  the  aid  of  Eqs.  (30)-(32)^ 
can  be  written  as 


where 


L*  ^  =  — V-  .  -  d.  .  .  -  e,  .  t,.  +  P, 

As2  1>J  i>J  i>:i  i>j  i’j 


(33) 


As2 


-li i± 

2As 


ft*  .  -  £*  ,  .]  , 

l  i+l,J  i-l, 3 j 


(34) 


and 


=  -A  .  ^  .  b  . 

1,3  Ai,j  *i,j  +  ®i,j  Zi,j 


+  tsS1  (**' -  ' 


(35) 
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where 


+ 


JuJ 

2As 


4W 


(36) 


=  fill 

2As 


(5 


i-2,j 


(37) 


The  number  of  unknowns  at  (i»j)  is  six,  namely:  (J>,  2^,  L^,  £,  2,  and 
L.S  Therefore,  Eqs.  (33)  and  (35)  must  be  supplemented  by  four  spline 
relations  to  close  the  system.  To  produce  a  compact  fourth-order  accurate 
system, s1  (4,0)  is  used  four  times:  (1)  to  relate  i p  to  2^,  (2)  to  relate 

to  L'f',  (3)  to  relate  £  to  2^,  and  (4)  to  relate  2*»  to  L** .  The  S^(4,0) 

tridiagonal  relationship  may  be  written  [12] 


AA, 


F.  .  .  + 


2  F 

0  h,i-i + 


BB 


Fi,j  + 


2  „F 


(1  +  o)  2 


i.j 


+  CCj  Fi,j+1  +  £i,j+l  "  °  ’ 


where  F  denotes  ip  or  or  (  or  and 


(38) 


,  _2o,  (2-W) 


J  An^d+O)  ’ 


(39) 


BB.**  2(1-0)  (1+0) ‘ 


(40) 


CC 


2(1  +20) 

j  ~  An.  x (1+0)0 


(41) 


An, 


o»o 


j  Anj-1  ’ 


(42) 


Anj  "  “j+i  “  "j 


(43) 
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The  number  of  unknowns  at  node  (i,i)  can  be  reduced  from  six  to  four 
by  substitution  of  the  expression  for  1$^  j,  Eq.  (33),  and  L^^j,  Eq.  (35), 

into  the  splines  relating  to  and  l ^  to  L^.  A  nonlinear  block  tridiagonal 
system  results  which  is  then  linearized  for  solution  by  the  Newton-Raphson 
technique.  If  superscript  n  denotes  the  interior  iteration  number,  associated 
with  the  nonlinear  spline-finite  difference  (SFD)  equations  on  line  i,  and  F 
denotes  one  of  the  four  unknowns  i|»,  then 

p(n+D  „  F(n)  +  6F(n)  ^  (44) 

where  6F^n^  is  assumed  to  be  a  small  correction.  The  linearized  system  is 
obtained  by  substituting  Eq.  (44)  into  the  nonlinear  system,  neglecting 

squares  of  6F^  and  writing  the  result  in  terms  of  the  corrections.  By 
solving  for  the  corrections  rather  than  the  functional  values  themselves, 
the  cancellation  errors  are  reduced  [14].  The  linearized  set  of  SFD  equations 
is  as  follows:  (the  n  superscript  is  understood.) 


(45) 
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2t’1’12'1-  S*l,l-1+  <“3  -  A.j-l’  <3-1  ■  <"2*l,i-l)  {C1.3-1 

As 


2b  ,(1+0) 2  r-  o  ~\  \b  F  2  1  r 

+  '[s2  — [  +  [BBj  difJJ  6£L  "  ei,jJ 


2b 

+  w  fi*i(J+i +  <ccj-di,-i+i)  6£t,j+r(ei,j+i)  6^i,j+i 

As 


D2,  . 
i,J  » 


2 


AA.  6*.^  +  +  BB.  6*ij+  (Ho)  «titJ 


+  CCj  +  64i,J+l  "  D3iJ 


2  xat 


liS?M-l  +  0“l.3-l  +  BBJ5'l.J  +  <14,,)  “m 


+  ccj  5cl,j+l  +  “i,J+l  ‘  °‘l.3  ’ 

The  right-hand  sides  of  Eqs.  (45) -(48)  are  given  by 

D1i,j'  A.j-1  -  (AAi  ■  °2ai.j-i>  *li-i 

3f  a2 

-  12Al—  lA  ~  tS)t ^  -  (1-W)2  *±fi  -  (BBj  -  a+C)\ti]  i\ 

-  ^  -  ^>ifJ  -  wUi  -  (“J  -  W  ll*« 


ah  -  *So 


t.j+i  . 


(49) 


-19- 


25  January  1982 
GHH:mmj 


and 


2  2  -i-i0 

D2  .  *  -  a  P.  .  .  -  (l+a)  P,  ,  -  P,  -  —  ii>.  .  , 

i.J  i,J-l  i,j  i.J+1  As2  ViJ-l 


-  (AA4  - 


o2d 


+  a2eJ 


i.j-1  i,J-l  l.J-l  i>j~l 


2b.  ,  (1+3)' 

-  ±jJ _ 

As2 


r±,i 


[BB^-  (1+0)' 


di,i] 


+ 

i,j 


(1+0) 2e 


i.j 


2b.  ... 

_ Li±l 

As2 


4> 


i,j+l 


'  <CCj  -  di,j+l> 


i,j+l  +  i,  j+1  4i,j+l 


and 


(50) 


D3 


.  =  -  AA.i|i.  .  .  -  O2^  .  .  -  BB.  <p.  .  -  (1-KJ) 2  if  . 
i.j  J  i.j-1  i .3-1  3  i,J  1,3 


~  ^i.J+l  "  ^i, j+1  ’ 


(51) 


and 


D4 


l.J  ‘  '  “j  C1,J-1  ‘  °  ‘l.J-l  "  BBJ 


C1,J  *  11+0,2  ‘l.J 


“  CCj  ?i,j+l  "  ^i, j+1 


(52) 


The  linearized  block  tridiagonal  system  for  the  corrections  may  be 
written  as  the  following  matrix  equation  along  line  i. 


-  A. 


Zj-1  +  Bj  Zj 


-  C 


j  j+1 


-  D. 


2  <  j  <  N 


(53) 


where  N  is  the  number  of  intervals  in  the  n  direction,  zj  is  the  four 
component  correction  vector  defined  by 


*1*,  6JlC]. 

1  *  J 


(54) 
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and  A,  B,  and  C  are  4x4  matrices  whose  elements  are  obtained  from  Eqs.  (45)- 
(48).  Finally,  Dj  is  a  four-component  column  vector  of  known  quantities 
with  components  (Dl,  D2,  D3,  D4)^  j  given  by  Eqs.  (49)-(52) . 

Spline-Finite  Difference  Equations  at  the  Wall 

The  boundary  conditions  on  the  wall,  Eqs.  (15a)  and  (15b),  in  SFD  form 


>J>,  ,  =  0  ,  (55 

1  =  0  .  (56 

Two  spline  relations  are  needed  to  close  the  system.  These  are  obtained  from 
Eq.  (15)  of  Rubin  and  Khosla  [12].  This  spline  relation,  obtained  from 
differentiation  of  S^(4,0),  is  a  two-point  formula  with  fourth-order 
truncation  error.  Applied  to  the  present  case  one  has  the  following  two 
equations : 

A  »1,2  -  *1.1>  -  ‘<1  +  <2>  -  Lt,l  *  0  •  <57 


4n  2  <Ci,2  "  '  Anl  +  "  Li,l  "  0  <5! 

The  spline  second  derivatives  L%  ^  aa<i  Li.l  are  eliminated  by  using  Eqs.  (33) 
and  (35)  evaluated  at  j  =  1. 

Substitution  of  Eq.  (35)  into  Eq.  (58)  gives  a  nonlinear  equation  which 
is  linearized  as  described  previously.  The  system  of  four  linear  equations 
at  the  wall  written  in  terms  of  corrections  is  as  follows: 

5*i.i  ■  D1i,i  •  (5 


Mf,i  *  B2i,i  • 


6  .  2bi,l 


2^-ff  «*i,l-  ArvT 

As  1 


— —  —  d.  .  +  e,  ,  6^.  ,  +  — — r  « 

\  i*1]  i>1  i*1  i*1  Anx2  i’2 


~  Wj  2  =  D3i  1 


(61) 
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and 


r  3f 

2As  ""1,1 


6i}/ 


1,1 


-  3f . 

j  +  - — C 
i,l  2As  t’i,l 


6  JL 


1,1 


An. 


3f .  ,  , 

-  +  i’l  J 

2  2As  ^i.l 


6C 


1,1 


-  3f .  .  . 

1,1  .  4 

Ai,l  +  2As  '♦'i.l  "  An. 


l 


6£*f  +  6?.  ,  -  ~~  6£*T  „  =  D4,  . 

i,!  A^2  i,2  Anx  1,2  1,1 


(62) 


The  right-hand  sides  of  Eqs.  (59)  and  (62)  are 


Dli,l  =  0  * 
D2i,l  =  0  • 


(63) 

(64) 


D3i,l  ’  Pi,l  + 


r  2b.  . 

6  +  i,l 


2  2 
An^  As4 


*i.i  + 


4"l  -  4t,l 


l.  ,  -  e.  .  C.  . 
i,l  i,l  i,l 


An/  ^.2  ^i  1,2  * 


(65) 


D4  «  B,  .  $  .  +  t,.  n  + 

i,i  l,i  1,1  £n  ^  i,l 


AnL  ~  Ai,l 


"i,! 


3f . 


+  2As  ^  "  ^i.l  “  2  Ci,2  +  An  *1,2 

An^  1 

The  quantities  Dli.,1  and  D2,  ,  are  both  zero  provided  the  correct  boundary 
conditions,  Eqs.  (55)  and  (56),  are  imposed  on  the  first  iterate  (n  *  1) . 


(66) 


The  linear  system,  Eqs.  (59) —(62) ,  may  be  written  as  the  following 
matrix  equation. 


B1  Z1  ~  C1  z2  "  D1 


(67) 
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where,  as  in  the  general  case,  and  are  4x4  matrices  with  known  elements 
obtained  from  Eqs.  (59)-(62)  and  is  a  four-component  column  vector  with 
components  given  by  Eqs.  (63)-(66). 

Spline-Finite  Difference  Equations  at  the  Outer  Edge 

The  boundary  conditions  on  the  outer  edge  of  the  computational  domain, 
Eqs.  (16a)  and  (16b),  in  SFD  form  are: 


^i.N+l 

"  *e.  * 

1 

(68) 

Ci,N+l 

=  «e.  • 

1 

(69) 

Two  spline  relations  are  also  needed  here  to  close  the  system.  Equation  (14)  of 
Rubin  and  Khosla  is  the  appropriate  relation,  analogous  to  the  spline  used  at  the 
wall .  Thus , 


and 


rfr  <*i,N+i  -  *i,N> 

N 


fi1?  +  2  tj 

AnN  l ,N  i,  N+lJ 


+ 


i,N+l 


=  0 


(70) 


pi,  N+l  "  Ci,N 


(71) 


Following  the  procedure  used  at  the  wall,  ,,,,  and  are  eliminated 

l,  N+l  l  ,N+1 

by  using  Eqs.  (33)  and  (35).  Then,  Eq.  (71)  is  linearized  and  all  four  equations, 
Eqs.  (68)-(71),  put  into  the  following  correction  form: 


and 


6^i,N+l  ~  D1i,N+l 


(72) 


and 


6^i,N+l  =  D2i,N+l 


(73) 


6  ~  „  + 


AnN2  ri*N  *’N 


2b 


i,N+l 


inN 


As 


6ip 


i,N+l 


f 


1%  +  di,N+l 


6X,i,N+l  "  ei,N+l  6;i,N+l  “  D3i,N+l 


(74) 


-23- 


25  January  1982 
GHHrmmj 


^aemm 


-24- 


25  January  1982 
GHH:mmj 


Then,  the  matrix  form  of  Eqs.  (72)-(75)  is 

"  \+l  ZN  +  BN+1  zn+l  =  °N+1  (80) 

where,  as  before,  and  BN+1  are  4x4  matrices 

with  known  elements,  from  Eqs.  (72)-(75),  and  DN+^  is  a  four-component  column 
vector  of  known  quantities  given  by  Eqs.  (76)-(79). 

Solution  of  Matrix  Equations 

Equations  (53) ,  (67)  ,  and  (80)  form  a  block  tridiagonal  system  which  may 
be  solved  by  standard  lower-upper  (L-U)  factorization  methods,  as  described  by 
Keller  (14].  To  avoid  buildup  of  numerical  errors  when  solving  the  subsidiary 
4x4  block  systems,  pivoting  must  be  used.  Blottner  [15]  has  developed  a  block 
tridiagonal  solver  which  performs  pivoting  within  the  blocks  and  also  minimizes 
storage  requirements  associated  with  the  L-U  decomposition  procedure.  The 
explanation  and  computer  code  for  this  solver  appears  on  pp.  27-32  of  Refer¬ 
ence  15.  Blottner' s  algorithm  has  been  adapted  to  the  solution  of  the 
present  system. 

Spline-Finite  Defference  Equations  at  a  Hap  Junction 

The  use  of  sheared  coordinates,  Eqs.  (7)  and  (8),  introduces  discontinuities 
into  the  metric  coefficients  along  an  s  =  constant  line  whenever  yw'(x)  and/or 
yw"(x)  are  discontinuous.  In  the  present  case,  only  yw"(x)  is  discontinuous 
and,  hence,  the  net  effect  on  the  SFD  equations  is  weak.  The  only  equation 
affected  is  the  stream  function  equation  because  central  differences  are  usei 
to  approximate  the  s-derivatives  which  at  a  map  junction  straddle  the  discon¬ 
tinuity.  Consequently,  the  SFD  approximation  of  the  stream  function  equation 
takes  on  a  special  form  at  map  discontinuities.  The  parabolized  vorticity 
equation  is  not  affected  because  it  contains  only  first  derivatives  in  s 
which  are  approximated  by  backward  differences  that  do  not  straddle  the 
discontinuity. 

The  map  junction  will  be  treated  in  a  general  manner  by  denoting  upstream 
and  downstream  regions  bordering  the  junction  by  superscript  (k)  and  (k+1) , 
respectively.  At  the  junction  smoothness  of  the  solution  requires  that 
and  its  various  derivatives  with  respect  to  x  and  y  be  continuous.  Note  that 
conditions  on  5  are  not  needed.  The  required  smoothness  conditions  are: 


^00  „  ^0*1)  .  , 

(81) 

<Vi<k)  -  <Vi<k+1>  ' 

(82) 

(«l>  )i(k)  “  (<p  )t(k+1)  , 

yy  1  yy  1 

(83) 

<Vi(k)  -  <Vi<krt)  • 

(84) 
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Using  the  chain  rule,  Eqs.  ( 81) - (84)  become 


i  )i(k)  (*  )± 

y  i  n  i 


2,(k) 


(k)  (‘^  (■=«) 

i  l  y  i  y  i 

(85) 

(<|)  ).(k)  =  Un  )  2](k+1)  OP  )  <k+1) 
nn  i  y  1  nni 

7 

(86) 

ip  ) . (k)  =  (s  i p  +  n  i p  )  (k+1) 
rn  i  x  rs  x  Tn  i 

(87) 

From  the  Appendix,  the  metric  coefficients  in  Eqs.  ( 85— ( 87)  are 


/I  a.  <2\1/2 

s  =  (1  +  y  ) 
x  w 


n  =  y  '  (y  -  y  )  (n  -  1)  , 

x  w  e  w 


n  =  (y  -  y  ) 
y  -'e 


At  the  map  junction,  y w  =  y w  and  y^'  =  yw'  ~  0.  Hence, 

i  i  i  i 

the  above  metric  cooef f icients  reduce  to 

.  <k>  - .  (k+1>  - 1  , 


„  <w  -  „  (k+1>  -  0  , 


X.  .  X.  . 

i.J  i.J 


„  00  _  n  (k+1) 

n  .  =  n  . 

yi  yi 


n  .  =  (y  -  y  ) 
yi  Je  wi 


As  a  result,  smoothness  relations  Eqs.  (85) — (87)  simplify  to 


0/>n)(k)  *  (^n)(k+1)  =  . 
n  i  n  i  11  i 


(\p  )(k)  =  0J|  )(k+1)  .  ) 

vynnyi  Ynn  1  rnn  i 


(90) 
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In  terms  of  spline  variables,  Eqs .  (88)  and  (89)  may  be  written 


£t(k)  -  £^(k+1)  =  t*  , 

m  i  -J  7 


L«k)  .  t*(W-l)  .  L*  .  (< 

i  i  i 

To  write  Eq.  (90)  in  discretized  form,  region  (k)  is  extended  one 
step  into  region  (k+1)  and  vice  versa  resulting  in  fictitious  node  lines. 
Unknowns  on  these  fictitious  node  lines  are  denoted  by  an  asterisk  (*) 
superscript.  Using  CD  approximations  for  in  Eq.  (90),  multiplying 
by  2As  and  rearranging  gives  the  following  result: 

^*(k+l)  +  *  |h  +  ^  (c 

^i-l,j  +  Vi+l,j  Vl.j  +  Vi+l,j  ’  U 

where  the  k  superscripts  have  been  omitted  from  values  of  tj>  lying  within 
their  respective  regions. 

The  values  of  on  the  fictitious  node  lines  appearing  on  the  left-hand 
side  of  Eq.  (93)  can  be  eliminated  by  using  the  stream  function  equation 
evaluated  at  (i,j),  the  node  point  common  to  both  regions  (k)  and  (k+1). 

Since  y  1  *  0  at  the  map  junction,  the  stream  function  equation  simplifies 


Ly  +  b  f  +  d  T  +  e  5  5  0 
rss 

Upon  using  a  CD  approximation  for  *j>ss,  Eq.  (94)  becomes 


*’3  As2  As2  ^  As2 


+  ,  $  X  +  e-  X  S,  X  -  0 

Let  the  i  index  at  the  junction  be  denoted  by  IJ  and  write  Eq.  (95)  at  IJ 
in  region  (k)  to  yield 

^U+l.j 

and  at  IJ  in  region  (k+1)  to  yield 


.*<*«> 

U-l.j 
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Substitution  of  these  results  into  Eq.  (93)  gives  the  following  special 
form  for  L1^  on  the  junction  line: 


I*  .  =  ]|>TT  .  -  d  £*  .  -  e  +  P 

IJ>3  Ag2  IJ.J  IJ.J  IJ,J 


where 


7  .  i  H  <k)  +  n<k+1) 

IJ.J  2  IJ.j  IJ,j 


(96) 

(97) 


(ij'x  T  i  •  +  ^ 

2  'rIJ-l,j  TIJ+l,j 


Equation  (96)  is  seen  to  be  the  same  as  the  general  case,  Eq.  (33),  provided 
d  is  replaced  by  d  and  P  is  simplified  as  above.  The  quantity  d  (related  to 
n^  and  hence  yw")  is  the  average  of  d  across  the  junction  line.  If  yw'  had 
been  discontinuous,  the  same  procedure  for  determining  the  special  expression 
for  L7  at  the  junction  would  have  been  applicable  but  a  more  complicated 
expression  would  have  resulted.  The  foregoing  treatment  of  junction  lines  is 
a  generalization  of  the  method  developed  by  Chmielewski  and  Hoffman  [16]. 

At  the  first  s  =  constant  line  downstream  of  a  map  junction,  the  second- 
order  backward  difference  used  for  s-derivatives  in  the  vorticity  equation 
requires  knowing  the  vorticity  on  the  fictitious  node  line  upstream  of  the 
junction.  To  avoid  the  problematical  calculation  of  this  value  of  vorticity, 
a  first-order  backward  difference  is  used.  Once  this  s-line  is  passed,  the 
scheme  reverts  to  using  second-order  backward  differences. 


Starting  Solution 


As  shown  in  the  schematic  in  Fig.  1,  the  boundary  conditions  on  the 
initial  line,  s  =  Sq,  (the  starting  line)  are  divided  into  two  zones: 

1.  Blasius  solution,  0  <_  y  £  6^ 

2.  Uniform  flow,  6q  <_  y  yg 

The  Blasius  solution,  which  is  assumed  to  hold  between  the  leading  edge 
(x  *  0)  and  the  initial  line  (x  *  xg) ,  is  thus  used  to  obtain  the  values  of 
ip,  l 5  and  ^  on  x  >  x~  for  0  <_  y  <_  6q.  The  Blasius  similarity  variables 
f  and  n  are  related  to  physical  variables  by 


(99) 
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where 


n 


(W!1/2 


N 


y 


The  other  three  required  quantities  are: 


(100) 


(101) 


C 


=  N 

N 


1/2 


f" 


(102) 


-  U  ff”  ,  (103) 

where  the  prime  denotes  differentiation  with  respect  to  n* 

The  Blasius  solution  is  calculated  as  follows:  The  similarity  form 
of  the  semi-infinite  flat  plate  boundary-layer  equation  is 

f",  +  ff«  =  o  f  (104) 

with  boundary  conditions, 

f=f'  =0onn=0  ,  (105) 

f'  =  i  on  n  =  .  (106) 

The  differential  equation  may  be  written  as  the  following  pair  of  differential 
equations : 


u  -  f'  ,  (107) 

u"  +  fu’  -  0  .  (108) 


Then,  upon  introduction  of  the  spline  variables: 
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the  pair  of  differential  equations  becomes,  at  point  j: 

Tj  ■  “j  •  (no) 

Lj  =  -  f  l  .  (Ill) 

Since  there  are  five  unknowns  at  j  (f,  £,  u,  £  and  L)  three  spline  relations 
are  needed  to  close  the  system.  To  be  consistent  with  the  parabolized 
solution,  sl(4,0)  is  used  three  times  and  the  number  of  unknowns  is  reduced 
to  f,  u  and  £  by  use  of  Eqs.  (110)  and  (111). 


In  the  uniform  flow  region  outside  the  boundary  layer,  the  initial 
flow  variables  are 


<J>  =  +  y  - 

0 


(112) 

(113) 

(114) 

(115) 


where  ^  is  the  value  of  the  stream  function  at  y  =  6q,  the  edge  of  the 
Blasius  boundary  layer. 


To  eliminate  the  requirement  of  two  initial  lines  to  start  the  solution, 
as  required  by  the  three-point  backward  difference  used  for  s-derivatives  in 
the  vorticity  equation,  a  two-point  backward  difference  is  used  at  the  first 
s  *  constant  line  downstream  of  s  =  Sg.  Thus,  the  solution  at  this  line  is 
only  first-order  accurate  in  As. 


Stream  Function  Equation  at  the  Downstream  Boundary 


In  region  II  the  coefficients  a,  c 


are  zero  because 


yw 


w 


=  0.  Hence, 


equation  ,  Eq.  22),  reduces  to 


L*  «  - 


b<|>  -  e  ? 


ss 


and  d,  given  by  Eqs.  (24),  (26)  and  (27), 
the  spline  form  of  the  stream  function 

(116) 


The  downstream  boundary 
was  (J>gs  ■  0  on  s  ■  s^. 
(i  =  I MAX)  is 


condition  on  that  terminates  the  elliptic  region 
Hence,  the  SFD  form  of  Eq.  (116)  on  this  boundary 


i y 

IMAX.j 


6 1  MAX  ^IMAXJ 


(117) 


-31- 


25  January  1982 
GHH:mmj 


From  the  Appendix,  Eq.  (12), 


(n  )  =  —  4>  y  *  i 

x  w  T  ■'w  ’ 

and,  hence,  using  Eq.  (13)  of  the  Appendix, 
2 


<v 


1  + 


4> 


~=l  +  v  '^  =  a^ 


From  the  Appendix,  Eq.  (9), 


(s  )  *  s  *  a 

x  w  x 


Then,  Eq.  (112)  may  be  written,  using  spline  notation,  as 

dC. 


dcp 
Re 

o  =  y, 


-  .  -  ,  ’  - - 

2  ds.  •'w  ds  w 


Equation  (123)  is  now  integrated  between  stations  i-1  and  i  to  yield 

i  i 

f  ' 

r  ^  1  I  r  d^.  ^ 


Re 

2 


c 

t  PW) 


wj 


i-1 


dC  'i 

y  ’  -“—I 

yw  ds 


ds 


i-1 


4  a  r  ds 
r  w 


i-1 


The  numerical  evaluation  of  d^^ds  can  be  avoided  by  integrating  the  first 
integral  on  the  right-hand  side  by  parts,  thus: 


r  <&„ )  i 

’  — *-l  ds  =  ly  ’  5  r 

w  i-1 


[yw  ds 
i-1 


• 

[y  " 

~  c 

a  w 

. 

ds 


i-1 


Then  the  pressure  coefficient  relation  becomes 


*p  p  Re 

rw  rw 

i  i-1 


(V  Cw}i 


1 

-  V  VH  -  \  o 

i-1 


ds 


A 


(123) 


(124) 
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where 


«  '  IT  +  *  0  *5  '  <125> 

The  final  result,  obtained  by  evaluating  the  integral  in  Eq.  (124)  by  the 
trapezoidal  rule,  is 


.  _L 

<V  Vi  -  <V  Vi-i  - 

[Qi-1 +  V 

(8i  -  Si-1> 

+  Re 

2 

i-1 

L 

(126) 


111.  RESULTS  AND  DISCUSSION 
Problem  Configurations 

Numerical  solutions  are  presented  for  six  cases  made  up  of  five  plate 
spacings  and  two  Reynolds  numbers  (10^  and  10^) .  The  parameters  for  each 
case  are  listed  in  Table  1.  In  cases  1-5,  the  rear  plate  is  positioned 
below  the  front  plate  so  that  an  adverse  pressure  gradient  occurs  on  the 
connecting  segment.  In  case  6,  the  situation  is  reversed  and  hence  a 
favorable  pressure  gradient  is  produced.  The  choice  of  the  number  of  grid 
lines  in  the  streamwise  (s)  direction  was  guided  by  the  work  of  Murphy  [9], 
who  solved  a  somewhat  similar  problem,  and  by  the  desire  to  keep  the  number 
of  mesh  points  at  a  minimum  so  that  run  times  would  not  be  excessive.  In 
Table  1,  the  guantities  NXj,  NXjj,  AND  NXjjj  refer  to  the  number  of  mesh 
intervals  in  regions  I,  II,  and  III,  respectively. 

In  the  present  problem,  Reynolds  number  has  a  somewhat  arbitrary 
meaning  because  by  rescaling  the  x  and  y  coordinates,  the  Reynolds  number 
is  correspondingly  rescaled  as  is  the  flow  field  since  it  is  laminar.  The 
important  point  to  keep  in  mind  is  that  the  viscous  layer  must  be  properly 
resolved  based  on  the  boundary  layer  scaling. 

Initial  Profile  Mesh  Choice 


As  already  mentioned,  the  initial  conditions  at  x  ■>  Xq  consist  of  a 
Blasius  boundary-layer  profile  adjacent  to  the  wall  and  an  outer  uniform 
profile.  A  rapidly  growing  mesh  in  y  is  used  which  is  small  near  the  wall 
to  properly  resolve  the  boundary  layer  and  is  large  in  the  outer  region 
where  flow  gradients  are  negligible.  This  vertical  grid  point  distribution, 
upon  conversion  from  y  to  n,  is  then  used  for  the  entire  computation  region. 
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Table  1 

Problem  Parameters 


CASE 

1A 

IB 

2 

3 

w 

5 

6 

Re 

io4 

io4 

io4 

io4 

io5 

io5 

io4 

xo 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

m 

0.20129 

0.20133 

0.20059 

0.20033 

0.20033 

0.20015 

0.30116 

x2 

0.30129 

0.30133 

0.30059 

0.30033 

0.30033 

0.30015 

0.40116 

xf 

0.40259 

0.40265 

0.40119 

0.40067 

0.60231 

yw 

1 

0.10 

0.10 

0.10 

0.10 

0.10 

0.10 

0.10 

yw 

W2 

0.085 

0.085 

0.090 

0.0925 

0.0925 

0.095 

0.11 

ye 

0.20 

0.20 

0.20 

0.20 

0.125 

0.125 

0.20 

NXX 

5 

10 

10 

10 

5 

5 

10 

NXn 

5 

10 

10 

10 

5 

5 

5 

NXI1I 

5 

10 

10 

10 

5 

5 

10 

As 

0.02026 

0.01013 

0.01006 

0.01003 

0.02007 

0.02003 

0.02012 
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The  nodal  distributions  in  zones  1  and  2  were  generated  by  geometric 
progressions  with  different  values  of  a  in  each.  In  zone  1  the  distribution 
of  the  Blasius  similarity  variable  rj  was  computed  based  on  Hoo  »  6  and  a  -  1.1 
which  gave  ten  intervals  across  the  boundary  layer.  The  present  eleven-point 
spline  solut.on  gives  f"(0)  =  0.46961  compared  to  f"(0)  =  0.46960  given  in 
Table  V.l  of  Rosenhead  [17].  From  rw  the  distribution  of  n  at  x  =  x„  in 
zone  1  is  computed  using  the  relation  J 


n 

j 


where 


(127) 


yj  = 


2x. 


1 1/2 


Re 


r|j  ,  for  0  <  ri  <_  6 


(128) 


In  zone  2  the  value  of  a  was  computed  because  yg  -  6q  and  the  number 
of  intervals  (taken  to  be  the  same  as  in  zone  1)  were  given.  Thus,  21  points 
spanned  the  computational  zone.  A  summary  of  the  important  initial  profile 
grid  parameters  appears  in  Table  2  for  the  two  Reynolds  numbers  used  in  this 
study.  As  can  be  seen,  the  values  of  and  a 2  are  nearly  the  same.  When 
02  was  taken  to  be  significantly  different  from  Oj_,  it  was  found  that  waviness 
appeared  in  the  flow  variables  in  the  essentially  inviscid  outer  region  and  no 
converged  solution  was  obtained.  The  quantity  n^  in  Table  1  is  the  value  of 

n  at  the  edge  of  the  initial  Blasius  profile  and  indicates  its  thickness 
compared  to  the  total  computational  zone  thickness  of  unity. 

Numerical  Solution  Procedure 


The  SLOR  procedure  was  initialized  by  writing  the  starting  profile,  as 
a  function  of  n,  at  all  s-stations.  Relaxation  sweeps  were  performed  in  the 
positive  s  direction  beginning  at  s  =  Sq  with  two  interior  iterations  of  the 
Newton-Raphson  procedure  at  each  s-station.  More  interior  iterations  per 
profile  were  tried,  but  convergence  of  the  entire  flow  field  was  not  improved. 
In  most  cases,  overrelaxation  was  performed  on  only  tp  and  JZT,  since  they  are 
princiaplly  governed  by  the  elliptic  stream  function  equation.  Above  a 
relaxation  factor  of  1.2  the  iterations  diverged.  Convergence  of  the  solution 
was  considered  attained  when 


max 


1 6z 


<  e 


where  z  denotes  the  four-component  column  vector  of  unknowns  defined  by  Eq.  (54) 
and  e  is  a  constant  0(1). 
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Numerical  Solutions 


The  numerical  solutions  presented  in  this  report  were  run  on  the  PSU 
IBM  370/3033  computer  using  double  precision  arithmetic.  The  two  flow 
parameters  presented  from  these  numerical  results  are  wall  friction 
coefficient  and  wall  pressure  coefficient.  The  wall  friction  coefficient, 
Cf,  is  related  to  the  wall  vorticity  by 


Cf  Re 

The  calculation  of  c  has  been  described  in  Section  II. 

Pw 

To  determine  the  effect  of  convergence  cutoff  on  the  numerical  solutions, 
case  1A  was  run  for  values  of  e  of  1.0,  0.1,  and  0.01.  The  results,  in  terms 
of  the  effect  on  C£,  are  presented  in  Table  3.  The  maximum  change  in  cf,  going 
from  e  of  1.0  to  0.01,  is  seen  to  be  at  most  three  in  the  fifth  decimal  place. 
Thus,  e  =  1.0  was  considered  to  be  satisfactory  for  all  cases.  A  cutoff  of 
unity  may  seem  unusually  large  until  one  realizes  that  the  convergence  test 
is  dominated  by  whose  maximum  magnitude  can  be  10^  or  10^. 

To  gain  a  rough  assessment  of  the  effect  of  the  horizontal  mesh  size 
on  the  solution,  two  values  of  As  (one  twice  the  other)  were  used  to  compute 
case  1  (  cases  1A  and  IB).  The  results  for  wall  friction  coefficient  are 
shown  in  Fig.  3.  The  differences  are  largest  in  the  unfavorable  pressure 
gradient  region  just  downstream  of  the  end  of  the  forward  plate  where  (Ac  ) 
is  about  0.002  at  x  =  0.24.  Other  than  the  comparison  of  the  initial  in'"‘x 

Blasius  solution  with  that  cf  Rosenhead,  no  assessment  has  been  made  of 
the  effect  on  the  overall  solution  of  changing  the  vertical  mesh  size  An. 

Overall  run  performance  for  the  six  cases  (and  sub-cases)  is  presented 
in  Table  4.  This  table  lists  number  of  node  points,  iterations  required 
for  convergence,  convergence  cutoff  and  CPU  time  required.  Cases  IB, 

2,  3,  and  6  had  nearly  twice  the  nodes  in  the  s-direction  as  the 
other  cases;  hence,  the  run  times  are  much  higher.  Run  time  can  undoubtedly 
be  reduced  by  careful  reprogramming  of  the  Blottner  solver  which  in  the 
present  application  recomputes  all  matrix  coefficients  at  every  mesh  point 
for  every  internal  (Newton-Raphson)  and  external  (relaxation)  iteration. 

Cases  1,  4,  and  5  all  contain  thin  separated  zones.  No  difficulties 
were  encountered  with  convergence  in  case  1.  To  obtain  convergence  at  the 
higher  Reynolds  number  (10^),  especially  in  case  4,  £  and  £.£  were  under¬ 
related.  The  relaxation  factors  used  were  0.5  in  case  4  and  0.7  in  case  5. 

The  relaxation  factor  for  ip  and  in  these  two  cases  was  1.0.  A  very  long 
separated  zone  occurred  in  case  4  which  extended  beyond  the  downstream 
boundary.  The  extent  of  this  separated  zone  is  reflected  in  the  larger 
number  of  iterations  required  for  convergence,  as  shown  in  Table  4  (71 
iterations  for  case  4  compared  to  26  iterations  for  case  5) .  The  stream¬ 
lines,  including  the  separation  streamline  and  the  extent  of  the  viscous 


1 
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Table  3 

Effect  of  Convergence  Cutoff  on  Wall  Friction  Coefficient,  Case  1A 


Sta 

— 

X 

e  =  1.0 

e  =  0.1 

e  =  0.01 

1 

0.10000 

.021002 

.021002 

.021002 

2 

0.12026 

.021174 

.021165 

.021164 

3 

0.14052 

.020821 

.020808 

.020807 

4 

0.16078 

.020545 

.020527 

.020526 

5 

0.18104 

.020842 

.020820 

.020817 

6 

0.20129 

.023059 

.023028 

.023025 

7 

0.22149 

.014322 

.014294 

.014291 

8 

0.24141 

.003817 

.003804 

.003803 

9 

0.26118 

-.000331 

-.000337 

-.000337 

10 

0.28110 

-.001570 

-.001573 

-.001573 

11 

0.30129 

-.001517 

-.001522 

-.001523 

12 

0.32155 

-.000222 

-.000233 

-.000234 

13 

0.34181 

.001290 

.001275 

.001274 

14 

0.36207 

.002787 

.002771 

.002769 

15 

0.38233 

.004161 

.004147 

.004145 

16 

0.40259 

.005414 

.005404 

.005403 

A 
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Run  Performance 
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Table  4 


Case  No . 

Nodes 

M  x  N 

Iterations 

To  Converge 

Cutoff 

(on 

CPU  * 

Time  (sec) 

1A 

16  x  21 

27 

1.0 

138 

1A 

16  x  21 

37 

0.1 

187 

1A 

16  x  21 

47 

0.01 

237 

IB 

31  x  21 

68 

1.0 

684 

2 

31  x  21 

53 

II 

531 

3 

31  x  21 

50 

M 

499 

4 

16  x  21 

71 

t» 

365 

5 

16  x  21 

26 

M 

132 

6 

26  x  21 

21 

M 

176 

*  IBM  370/3033 
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region  (determined  as  (  >  0.01),  are  plotted  in  Fig.  4  for  case  1  and 
Fig.  5  for  case  4.  These  plots  show  that  the  viscous  region  covers 
about  half  the  total  computational  zone. 

Numerically  calculated  distributions  for  wall  friction  coefficient 
for  all  cases  are  shown  in  Figs.  6-11  and  for  wall  pressure  coefficient 
for  cases  1,  2,  3,  and  6  in  Figs.  12-15.  At  the  lower  Reynolds  number 
(10^*)  comparisons  are  presented  in  Figs.  6-8  and  11-15  with  potential 
flow-boundary  layer  solutions  having  no  strong  interaction  coupling. 

These  comparisons  clearly  show  the  extent  and  magnitude  of  the  strong 
interaction  effect.  The  potential  flow-boundary  layer  solutions  were 
obtained  by  a  modification  of  the  SFD  parabolized  vorticity  code. 

In  cases  1,  2,  and  3,  where  the  rear  plate  is  below  the  front  plate, 
the  flow  first  accelerates  slightly  as  it  approaches  the  convex  center 
wall  section  and  then  abruptly  slows  down.  This  behavior  is  demonstrated 
by  the  wall  pressure  distribution,  as  shown  in  Figs.  12-14.  The  potential 
flow  is  seen  to  overpredict  the  levels  of  the  favorable  and  unfavorable 
pressure  gradients.  Consequently,  the  peak  in  Cf,  corresponding  to  the 
favorable  pressure  gradient,  is  much  too  high,  as  shown  in  Figs.  6-8. 

The  magnitude  of  the  potential  flow  adverse  pressure  gradient  is  so  large 
that  the  boundary  layer  separates  in  all  of  these  cases  even  though  the 
parabolized  vorticity  solution  separates  (and  reattaches)  only  in  case  1, 
as  shown  in  Fig.  6.  The  boundary-layer  separation  points  in  Figs.  6-8 
are  not  the  exact  separation  points,  but  the  station  where  the  boundary- 
layer  calculation  procedure  failed  to  converge  which  indicates  the  procedure 
has  marched  across  a  separation  point. 

In  case  6,  where  the  rear  plate  is  above  the  front  plate,  the  boundary- 
layer  solution  separated  just  upstream  of  the  beginning  of  the  curved  mid¬ 
section,  as  seen  in  the  plot  of  c^  in  Fig.  11.  The  reason  for  this  behavior 
is  shown  in  the  wall  pressure  distribution.  Fig.  15,  where  a  strong  adverse 
pressure  gradient  is  predicted  by  the  potential  flow  solution  in  that  region 
The  flow  is  slowing  down  as  it  encounters  the  concave  wall.  The  parabolized 
vorticity  solution  does  not  separate,  however. 

The  wall  pressure  distributions,  Figs.  12-15,  clearly  show  that  the 
large  viscous  region  reduces  the  peaks  and,  hence,  the  magnitude  of  the 
gradients.  Consequently,  only  in  case  1  is  the  adverse  pressure  gradient 
strong  enough  to  cause  separation.  With  the  strong  interaction  effects 
cutting  down  on  the  wall  pressure  peaks,  the  peaks  in  wall  friction 
coefficient  are  also  greatly  reduced. 

IV.  CONCLUSIONS 


For  the  present  test  geometry  and  Reynolds  numbers  the  parabolized 
vorticity  approximation,  using  a  spline-finite  difference  discretization,  is  an 
efficient  as  well  as  effective  way  of  treating  the  strong  interaction  problem, 
attached  flows,  convergence  of  the  SLOR  procedure  is  rapid.  The  procedure 


For 
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Is  also  able  to  predict  thin  separated  zones  with  acceptable  convergence 
properties . 

The  results  of  this  study  reinforce  the  well-known  observation  that 
to  obtain  meaningful  numerical  results  at  high  Reynolds  numbers  a  viscous 
layer  adjacent  to  a  wall  must  be  properly  resolved  according  to  the 
boundary-layer  scaling.  Proper  resolution  with  a  relatively  sparse  grid 
is  achieved  by  use  of  fourth-order  accurate  splines  to  approximate  flow 
quantities. 

The  effect  of  the  slightly  non-ortho gonal  coordinate  system  on  the 
accuracy  of  parabolized  vorticity  model  solutions  will  be  dealt  with  in 
a  future  report. 
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Fig.  3.  Effect  of  Streetwise  Scepsis.  on  W1  Friction  Coefficient,  Case  1,  Reynolds  No 
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Fig.  11.  Wall  Friction  Coefficient,  Case 
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Fig.  13.  Wall  Pressure  Coefficient,  Case 
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Fig.  14.  Wall  Pressure  Coefficient,  Case 
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APPENDIX:  METRIC  COEFFICIENTS 


The  sheared  coordinate  transformation  is 


s  =  ct(x)  dx  ,  (A.l) 


n  =  <j>(x)  •  y  -  A(x) 


(A. 2) 


where 


a 


(l  +  y 


.2  1/2 

w 


(A.  3) 


4> 


X 


(A.  A) 


(A.  5) 


Thus,  a,  <(>,  and  X  depend  on  the  wall  equation  y  (x)  and  its  derivative  y’(x). 
The  wall  equation  is  as  follows:  In  regions  I  and  III 


x0  ixi 


x2  <  x  <  xf 


(A6a) 


(A6.b) 


while  in  region  II,  y  (x)  is  given  by  a  cubic  such  that  yw  and  y  '  are 
continuous  at  the  end  points  x^  and  X£.  The  resulting  equation  is 


yw(x)  =  ywL  +  <w2  ‘  yw1)  "  25) 


(A.  7) 
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